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ABSTRACT 

The  multidimensional  inverse  scattering  problem  for  an  acoustic  medium  is  con¬ 
sidered  within  the  homogeneous  background  Born  approximation.  The  medium  is 
probed  by  wide-band  plane  wave  sources,  and  the  scattered  field  is  observed  along 
straight-line  receiver  arrays.  The  objective  is  to  reconstruct  simultaneously  the 
velocity  and  density  profiles  of  the  medium.  The  time  traces  observed  at  the  re¬ 
ceivers  ?ire  appropriately  filtered  to  obtain  generalized  projections  of  the  velocity 
and  density  scattering  potentials,  which  are  related  to  the  velocity  and  density  per¬ 
turbations  of  the  medium  with  respect  to  their  nominal  values.  The  generalized 
projections  are  weighted  integrals  of  the  scattering  potentials;  in  two  dimensions 
the  weighting  functions  are  concentrated  along  parabolas,  in  three  dimensions  they 
are  concentrated  over  circular  paraboloids.  The  reconstruction  problem  for  the 
generalized  projections  is  formulated  in  a  way  similar  to  the  problem  of  x-ray,  or 
straight-line  tomography.  The  solution  is  expressed  as  a  backprojection  operation 
followed  by  a  two  or  three  dimensional  space-invariant  filtering  operation.  In  the 
Fourier  domain,  the  resulting  image  is  a  linear  combination  of  the  velocity  and  den¬ 
sity  scattering  potentials,  where  the  coefficients  depend  on  the  angle  of  incidence  of 
the  probing  wave.  Therefore,  two  or  more  different  angles  of  incidence  are  necessary 
to  reconstruct  the  velocity  and  density  scattering  potentials  separately. 
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INTRODUCTION 

In  this  paper,  we  consider  the  multidimensional  inverse  scattering  problem  for 
an  acoustic  medium  within  the  homogeneous  background  Born  approximation.  An 
acoustic  medium  is  probed  by  wide-band  plane  wave  sources,  and  the  scattered 
field  is  observed  along  straight-line  receiver  arrays.  The  objective  is  to  reconstruct 
simultaneously  the  velocity  and  density  profiles  of  the  medium. 

The  1-D  velocity/density  profile  inversion  problem  has  been  studied  by  a  number 
of  researchers,  including  Raz  (1981a),  Coen  (1981),  Hooshyar  and  Razavy  (1983), 
Yagle  and  Levy  (1984).  This  problem  can  in  principle  be  solved  exactly. 

The  multidimensional  problem,  where  the  velocity  and  density  profiles  are  al¬ 
lowed  to  vary  in  two  or  three  dimensions,  has  also  interested  several  researchers. 
For  the  multidimensional  problem,  no  exact  solution  exists.  Several  approximate 
inversion  techniques  have  been  proposed,  which  linearize  the  scattering  integral 
equations  by  using  the  Bom  or  Rytov  approximations.  In  this  context,  it  was 
shown  that  the  experimental  requirements  of  single  parameter  inversion  problems, 
where  the  medium  density  is  constant  and  only  the  velocity  varies,  and  of  mul¬ 
tiparameter  inversion  problems,  where  both  the  density  and  the  velocity  need  to 
be  reconstructed,  are  different.  For  the  single  parameter  case,  a  single  scattering 
experiment  is  sufficient  to  -  at  least  partially  -  reconstruct  the  object  of  interest, 
whereas  several  experiments  are  necessary  for  multiparameter  problems. 

Raz  (1981b),  and  Clayton  and  Stolt  (1981)  have  solved  the  multidimensional 
problem  for  an  experimental  setup  where  sources  and  receivers  are  available  at  all 
points  on  a  plane,  and  scattered  waves  are  measured  at  all  frequencies.  Coen, 
Cheney  and  Weglein  (1984),  and  Ramm  and  Weglein  (1984)  have  solved  the  same 
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problem  for  a  complete  set  of  sources  and  receivers,  but  using  only  two  temporal 
frequencies.  Hooshyar  and  Weglein  (1986)  have  used  two  wide-band  point  sources, 
where  the  distance  between  the  sources  is  required  to  be  small  compared  to  the 
distance  from  the  sources  to  the  scatterers.  Beylkin  and  Burridge  (1987)  have 
presented  a  solution  for  a  medium  with  a  variable  background,  with  sources  and 
receivers  surrounding  the  medium. 

The  class  of  problems  where  the  medium  is  probed  by  plane  waves  has  been 
investigated  by  Norton  and  Devaney.  Norton  (1983)  has  used  a  flat  transducer  as  a 
source  of  wide-band,  plane  wave  illumination,  and  as  a  receiver  of  the  backscattered 
waves.  A  second  transducer,  oriented  at  a  different  angle  with  respect  to  the  first, 
is  used  as  a  receiver  only.  The  two  transducers  are  rotated  together  180°  around 
the  object,  and  the  scattered  waves  are  recorded  at  all  angles  during  the  rotation. 
Devaney  (1985)  has  extended  the  diffraction  tomography  theory  to  the  variable 
density  case.  In  this  work  the  transmitted  waves  are  measured  on  a  plane  parallel 
to  the  incident  plane  wave  front,  and  the  experiment  is  repeated  for  all  view  angles. 
In  this  scheme,  two  temporal  frequencies  are  used  in  the  insonifying  wave. 

The  present  paper  is  a  generalization  of  a  previous  work  of  the  authors  (Ozbek 
and  Levy,  1987),  where  the  multidimensional  inverse  scattering  problem  for  a  con¬ 
stant  density  acoustic  medium  was  formulated  and  solved  as  a  generalized  tomo¬ 
graphic  problem.  In  this  paper,  we  similarly  filter  the  time  traces  observed  at  the 
receivers  to  obtain  generalized  projections  of  the  velocity  and  density  scattering  po¬ 
tentials,  which  are  related  to  the  velocity  and  density  variations  in  the  medium.  The 
generalized  projections  are  weighted  integrals  of  the  scattering  functions;  in  the  two- 
dimensional  geometry  the  weighting  functions  are  concentrated  along  parabolas,  in 
the  three-dimensional  geometry  they  are  concentrated  over  circular  paraboloids. 
Thus  the  inverse  scattering  problem  is  again  posed  as  a  generalized  tomographic  or 
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integral  geometric  problem. 

The  reconstruction  problem  for  the  generalized  projections  is  formulated  in  a 
way  similar  to  the  problem  of  x-ray,  or  straight-line  tomography.  The  solution 
is  expressed  as  a  backprojection  operation  where  we  sum  the  contributions  of  all 
projections  passing  through  a  given  point  in  space,  followed  by  a  two  or  three 
dimensional  space-invariant  filtering  operation.  In  the  Fourier  domain,  the  resulting 
image  is  a  linear  combination  of  the  velocity  and  density  scattering  potentials,  where 
the  coefficients  depend  on  the  angle  of  incidence  of  the  probing  wave.  Therefore, 
two  or  more  different  angles  of  incidence  are  necessary  to  solve  for  the  velocity  and 
density  scattering  potentials  separately. 

The  main  difference  between  the  approach  that  we  propose  here  and  the  diffrac¬ 
tion  tomography  technique  developed  by  Devaney  (1985)  is  that  we  use  wide-band 
plane  waves  at  just  a  few  angles  of  incidence  to  reconstruct  the  medium  density 
and  velocity,  whereas  the  diffraction  tomography  formulation  relies  on  narrowband 
plane  wave  sources  at  just  two,  and  in  practice  several,  frequencies,  but  for  all 
angles  of  incidence.  These  two  approaches  are  in  some  sense  dual  of  each  other, 
since  they  trade  wavevector  orientation  against  frequency.  This  distinction  leads  to 
significantly  different  processing  algorithms,  and  in  fact,  as  mentioned  above,  the 
inversion  procedure  proposed  here  is  significantly  closer  to  x-ray  tomography  than 
to  diffraction  tomography. 

Also,  the  different  choice  of  experimental  conditions  appearing  here  is  a  reflection 
of  the  difference  existing  between  medical  imaging  applications,  which  motivated 
the  diffraction  tomography  approach,  and  geophysical  tomography  problems,  which 
are  at  the  origin  of  the  present  work.  In  medical  imaging  it  is  possible  to  rotate  the 
object  to  be  imaged,  i.e.  the  patient,  over  a  360°  range,  but  in  exploration  geophysics 
only  a  very  limited  range  of  angles  of  incidence  can  be  achieved  for  surface,  or  even 
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vertical  seismic  profiling  recordings. 

The  paper  is  organized  as  follows:  we  treat  the  2-D  case  in  detail  in  Sections 
I-III,  and  just  summarize  the  results  for  the  3-D  case  in  Section  IV.  In  Section  I, 
the  inverse  scattering  problem  is  formulated  within  the  Born  approximation  and 
redefined  as  a  generalized  tomographic  problem.  The  backprojection  operation  is 
defined  and  related  to  the  generalized  projections  in  Section  II.  In  Section  III,  the 
separate  reconstruction  of  the  velocity  and  density  scattering  potentials  is  discussed. 
In  this  context,  it  is  shown  that  substantially  more  than  two  plane-wave  experiments 
are  required  in  order  to  be  able  to  recover  the  velocity  and  density  perturbations 
in  a  numerically  reliable  way.  We  summarize  the  results  for  the  3-D  geometry  in 
Section  IV.  A  2-D  numerical  example  is  presented  in  Section  V,  and  Section  VI 
contains  some  conclusions. 

L  PROBLEM  DESCRIPTION 

In  this  paper  we  will  treat  the  two-dimensional  case  in  detail,  and  summarize  the 
results  for  the  three-dimensional  case.  Consider  the  scattering  experiment  described 
in  Fig.  1.  A  2-D  acoustic  medium  is  probed  by  a  wide-band  plane  wave  and  the 
scattered  field  is  observed  along  a  straight-line  receiver  array.  The  Fourier  transform 
of  the  pressure  field  at  point  x  =  {x,y)  satisfies  (Chernov,  1960) 

^ 

where  c(^  is  the  propagation  velocity,  and  /)(z)  is  the  density  of  the  medium  at 
point  X.  The  acoustic  equation  (1)  can  be  rewritten  as 


(V^  +  P)P{x,uj)  =  -k^U,{x)P{x,u})  +  •  VP{x,a;), 


(2) 
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where 

(3) 

Up{x)  =  In  ^  \ 

Po 

(4) 

and  k  =  w/co  is  the  wavenumber.  Here,  Ue{^  and  Uf{^  are  respectively  called 
the  velocity  and  the  density  scattering  potentials.  Cg  and  Po  are  respectively  the 
velocity  and  the  density  of  the  backgroimd  medium.  We  assume  that  c(i)  and  p(^ 
do  not  deviate  significantly  from  their  nominal  values  Cg  and  pg;  consequently  17c(x) 
and  17p(x)  are  small  with  respect  to  1.  We  also  assume  that  17c  ($)  and  17p(x)  have  a 
bounded  support  “V ,  which  is  located  completely  on  one  side  of  the  receiver  array. 

The  probing  wave  Pg{x,uj)  satisfies 

(V^  +  fc2)Pc(^u;)=0,  (5) 

so  that  the  scattered  field  P,(x,  u)  =  P(xj  u;)  —  Pg{^  u)  obeys 

(V*  +  k^)P,{x, w)  =  -k^Ug{x)P[x,  w)  +  V17p(x)  •  VP(x,  w).  (6) 

The  solution  of  (6)  is  given  by  the  Lippmann-Schwinger  equation  (Taylor,  1972) 

[k'^Ug{^)P{^,i^)  -  V17p(x')  .  VP(x',a;)]  C?c(x,x',a;),  (7) 

where  Go(x,x',a;)  is  the  Green’s  function  associated  with  a  point  source  in  a  homo¬ 
geneous  medium: 

(V*  +  A:^)Go(x,x',a;)  =  -<5(x- x').  (8) 

Equation  (7)  demonstrates  the  nonlinear  relation  that  exists  between  the  po¬ 
tentials  Ue{^,  Up{x)  and  the  pressure  field  P(^u)).  To  linearize  this  equation  we 
adopt  the  Born  approximation,  whereby  we  assume  Ps(x,cj)  <C  Po(^w).  Then  the 
Lippmann-Schwinger  equation  becomes 

PsU,c^)  «  / dx!  [k^l7g(x!)Pg(x',u:)  -  VU,(x')  ■  VPp(x',c*;)]  Gg(x,x',co).  (9) 
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The  Born  approximation  assumes  that  the  scattered  field  Pg  cj)  is  small  through¬ 
out  the  volume  of  the  object,  which  requires  that  both  the  magnitude  of  the  scat¬ 
tering  potentials  Ue  and  Up  be  small  and  that  the  volume  of  the  region  V  be  small 
with  respect  to  the  dominant  wavelength  of  the  probing  wave. 

If,  instead  of  the  Born  approximation  one  uses  the  Rytov  approximation,  the 
requirement  that  the  size  of  the  scattering  region  be  small  can  be  relaxed  (Chernov, 
1960;  Tatarski,  1961;  Devaney,  1981).  The  Rytov  approximation  is  obtained  by 
representing  the  total  pressure  field  P(^a;)  in  terms  of  its  complex  phase  and  lin¬ 
earizing  the  resulting  Riccati  equation  satisfied  by  the  phase  fluctuation  (Devaney, 
1985). 

On  the  other  hand  the  Bom  approximation  is  more  accurate  than  the  Rytov 
approximation  for  reflected  waves  (Beydoun  and  Tarantola,  1986).  For  the  setup 
considered  in  this  paper,  where  the  wideband  property  of  the  probing  wave  replaces 
the  large  number  of  view  angles  available  in  diffraction  tomography,  it  was  shown 
in  Ozbek  and  Levy  (1987)  that  the  reflected  wave  configuration  provides  the  most 
coverage  in  the  Fourier  domain  for  a  bandlimited  source.  Therefore,  we  adopt  the 
Born  approximation  in  this  paper,  although  similar  results  can  also  be  derived  for 
the  Rytov  approximation. 

Next,  we  simplify  the  second  term  in  the  integrand  in  equation  (9)  by  using  the 
identity  (Norton,  1987) 

(VUp  .  VP,)Gp  =  V  ■  {UpGpVPp)  -  UpV  •  {GpVPp),  (10) 

and  applying  the  divergence  theorem.  Over  a  surface  S  located  outside  the  domain 
T  where  the  density  inhomogeneities  are  located,  we  have  Up{^  —  \n{p{^/ po)  —  0, 
so  that 


/  dPV  >[Up{P)G,{x,^,w)VPo{P,uj)] 

j  y 
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=  dsUf, (2') Go (x,  s' ,  w) n {5)  •  V Po ( s' ,  w) 

=  0.  (11) 

For  the  second  term  in  equation  (10),  we  have 

V  .  (GoVPo)  =  VGo  •  VPo  +  GoV^Po  =  VGo  •  VPo  -  k'^G^Po-  (12) 

The  incident  wave  is  given  as 

Pa(x',a;)  =  e‘'=^-^',  (13) 

where  £  =  (cos0,  sin0)  is  the  unit  vector  which  indicates  the  angle  of  incidence 
of  the  plane-wave  source.  Therefore,  within  the  Born  approximation,  the  scattered 
field  at  a  receiver  point  ^  is  given  by 

Pt  (i,  w)  =  y  dP  [k^\Uc (s')  -Uf, (s') ] Po (s' ,  w) Go (1^,  s',  w) 
+i7,(s')V,.Po(s',w) .  V..Go(i,s',  w)} 

=  J  dp[k'^[Uo{P)  -  P’p(s')]Go(|,s',a;) 

+ikU,[P)[h  V,.Go(i,s',  w)]}  Po(s',a;).  (14) 

To  compare  the  far-field  scattering  patterns  that  are  due  to  velocity  and  density 
perturbations,  let  us  use  the  far  field  approximation  fe(s'  —  »  1.  We  find 

V*»Go(£,s',a;)  «  ik^, — ||-Go(^,s',  w),  (15) 

"  “  IS.  -  ii  “ 

which  is  valid  for  both  the  2-D  and  3-D  Green’s  functions.  Then  (14)  becomes 

Ps{i,<^)  «  k^  J  dx’[Uo{P)  -  [1  +  cosa(|,x',£)]?7p(s')}  Go($,x',cj)Po(x\u),  (16) 

where  a;(£,  s',  £)  is  the  angle  between  the  vectors  — £  and  ^  —  x'.  From  a  physical 
point  of  view,  o:(^,s',£)  is  the  angle  at  point  x'  between  the  ray  linking  s'  to  the 
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source,  which  is  specified  by  vector  since  the  incident  wave  is  a  plane  wave, 
and  the  ray  linking  z'  to  receiver  Since  the  background  medium  is  assumed  to 
be  constant,  these  rays  are  straight  lines.  The  different  weights  in  front  of 
and  in  (16)  correspond  to  different  scattering  patterns  which  are  plotted  in 

Figs.  2a  and  2b.  As  indicated  by  Fig.  2a,  the  scattering  pattern  due  to  Uc{x)  is 
that  of  a  monopole,  whereas  Fig.  2b  shows  that  the  scattering  pattern  due  to 
is  that  of  the  sum  of  a  monopole  and  a  dipole.  Therefore,  the  scattering  due  to 
density  perturbations  is  most  prominent  for  reflected  waves,  and  least  prominent 
for  transmitted  ones. 

Equation  (14)  can  be  written  in  a  slightly  different  form  if  we  introduce  the 
scattering  potential 

(17) 

Ko 

associated  to  the  compressibility  /c(z)  =  l/c^(^p(z).  We  find 

=  In  -  U,[x)  «  U,{x)  -  (18) 

c  \x) 

for  c(z)  near  Co,  as  assumed  above.  Therefore  (14)  can  be  written  as 

(19) 

and  the  objective  of  this  paper  will  be  to  solve  this  integral  equation  for  the  com¬ 
pressibility  and  the  density  scattering  potentials  C/k(z)  and  ?7p(z),  or  equivalently 
for  Up(^  and  the  velocity  scattering  potential  f7c(z). 

For  the  2-D  geometry  under  consideration,  the  Green’s  function  is  given  by 

=  ^hS,^^(A:|z  -  z'l),  (20) 

where  Ho^^(')  indicates  the  Hankel  function  of  order  zero  and  type  one.  In  the 
following,  it  will  be  assumed  that  the  receivers  are  located  along  a  straight  line 
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perpendicular  to  the  unit  vector  ^  =  (cos^,  sin<^)  and  whose  distance  from  the 

A 

origin  in  the  direction  <f>  is  p,  as  shown  in  Fig.  1.  The  position  of  an  arbitrary  receiver 

A  aJ.  aX 

along  this  line  is  therefore  given  by  £  =  ,  where  ^  =  (sin^,  —  cos^)  is  a 

unit  vector  perpendicular  to  and  ^  is  an  arbitrary  coordinate.  Then  (19)  can  be 
expressed  as 


where 


=  FU,k)  ^  F4^,k)  -  F,i^,k), 


t{^,k)  =  y/  -  i\], 

=  -  1/  dx'U,{x;)e''‘^-^‘ 


(21) 

(22) 


s'  - 


HS‘>(4|x'-il),  (23) 


and  where  =  — dH^*^{u)/dv  is  the  Hankel  function  of  order  1  and  type  1. 

We  define  the  inverse  Fourier  transform  of  F{^,k)  with  respect  to  k  as  g{C,r): 

aU: r)^Lr  dkP(i, (24) 

Z7C  J  — oo 

F)  and  gp{i,  r)  are  similarly  defined  as  the  inverse  Fourier  transforms  of  F^{^,  k) 
and  Fp{C,k),  respectively.  Taking  into  account  the  fact  that  (see  Morse  and  Fesh- 
bach,  1953,  pp.  1362-1363) 

l(r  —  u) 


■'{fHW(M}  = 


V~r 


where  l(-)  is  the  unit  step  function,  and  the  relationship 

that  follows  from  (25) ,  we  find  that 


(26) 


(26) 


9(£.’')  =9<(£,r)  -9,(f,r), 


(27) 
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where 


l(r  -  £  •  x'  -  -  ||) 


■'  V  (>*-£•  s' 


and 


- 1^'  -  ir 

r  —  0  •  s'1  l(r  —  £ •  s'  —  |s'  —  £|) 


(28) 


.'-iiJLk'-ii 


\l{r  -  £*x')2  -  lx'  -  £| 


(29) 


Equation  (28)  expresses  fif,c(^5^)  as  a  weighted  integral  of  the  compressibility 
scattering  potential  Uk(x)  where  the  weighting  function  is  nonzero  in  a  region  with 
parabolic  support.  The  parabola  satisfies  the  equation  r  =  £  •  s  —  |s  —  £|  where  r , 

A  A 

£  and  ^  are  given  and  s  varies.  The  directrix  of  this  parabola  is  the  line  £  •  x  =  r 

A 

which  is  perpendicular  to  the  direction  £  of  incidence  of  the  probing  wave,  and 
whose  focus  is  the  receiver  point  £.  The  weighting  function  becomes  infinite  for 
values  of  s  along  this  parabola,  so  that  the  largest  contribution  to  the  integral  is 
made  by  the  values  of  17,t(s)  which  lie  along  the  parabola.  In  some  sense,  gK(£,r) 
is  a  projection  of  with  respect  to  a  function  whose  singularities  are  algebraic 

and  located  along  a  parabola. 

In  (29),  gp(^,r)  is  similarly  expressed  as  a  weighted  integral  of  the  density  scat¬ 
tering  potential  U^{^.  The  weighting  function  again  has  a  parabolic  support,  but 
it  contains  two  additional  factors.  The  first  factor,  £  •  (s'  —  £)/|s'  —  ^|,  is  again 

A  A  A 

identified  as  cosq:(£,s',£),  where  q:(^,s',£)  is  the  angle  between  the  vectors  — £  and 
£—  s'.  The  second  factor  in  (29),  (r  —  £•  x')/|s'  —  £|,  is  a  term  which  equals  unity 
on  the  parabolic  wavefront  of  the  weighting  function,  and  grows  as  the  focus  of  the 
parabola  (the  receiver  point)  is  approached. 

In  the  following,  it  will  be  assumed  that  the  projections  g{^,r)  constitute  the 
data  obtained  from  a  single  plane-wave  scattering  experiment.  To  see  why  this  is 
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the  case,  observe  from  (21)  and  (24)  that 

Then  since 

^{1(0}  =  -^  +  7r(5(u;),  (31) 

lU) 

and  observing  from  (21)  that  P,(^,a;  =  0)  =  0,  we  find  that 

9i^,r)  =  -27rCo  I  f  ^  dr  f  dsP,{^,s)  -  ^  f  dr  j  dsP,(^,5)| .  (32) 

00  J —OO  Z  J --00  J  ^oo  J 

The  projections  g(^,  r)  are  therefore  obtained  by  integrating  twice  the  time  domain 
scattered  field  Fg(^,t)  observed  at  and  then  subtracting  a  constant  equal  to  half 
the  value  of  the  double  integral  at  f  =  +oo.  This  shows  that  the  knowledge  of  the 
observed  scattered  field  P,(^,  t)  for  all  t  is  equivalent  to  that  of  projections  g{^,  r) 
for  all  r. 

On  the  basis  of  the  above  observations,  the  inverse  scattering  problem  can  be 
formulated  as  follows:  given  the  generalized  projections 

{^(^j’')  :  ~oo  <  ^  <  OO,  —CO  <  r  <  oo}, 

we  want  to  reconstruct  the  scattering  potentials  17^ (^  and  Up(x). 


II.  BACKPROJECTION  OPERATION 


Proceeding  as  in  the  constant  density  inversion  problem  treated  in  Ozbek  and 
Levy  (1987),  the  first  step  of  our  inversion  procedure  is  to  perform  a  backprojection 
operation  on  the  projections  ff(^,r).  We  define  it  as 


rr  ,  \  A  J  fy.  \  l{r-9-X-  X-£) 

Ub{x)  =  dU  drg{i, r) -= . - . 

V(r-£'z)2  - 


(33) 
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At  a  given  point  x,  this  operation  sums  the  contributions  of  x  to  the  projections 
Sf(^,r),  as  appearing  in  equations  (27)-{29).  By  performing  this  backprojection 
operation  for  every  point  in  the  plane,  this  gives  an  image,  Ub{x).  Ub{^  is  the 
analog  of  the  backprojected  image  obtained  in  x-ray  tomography  (Deans,  1983)  by 
summing  the  line  projections  passing  through  a  given  point. 

Our  first  objective  is  to  relate  Ub{^  to  the  projections  g[i,T)  in  the  frequency 
domain.  It  can  be  shown  that  the  2-D  Fourier  transform  of  Ub{^  is  given  by  (see 
Ozbek  and  Levy,  1987,  Section  3  and  Appendix  A) 


UbHc)  =  I  dxUBi^e-'^-^ 


ITT 

in 


ki:  k=-^] 

2i-r  ^  2k’§)’ 


(34) 


where 

^  /‘OO  /*oo 

g{k^,kr)=  di  I  dr5(^,  (35) 

is  the  2-D  Fourier  transform  of  g{^,r),  k=  {kx,ky),  k  =  |fc|,  A  =  (fc^  -  kl,  2kxky), 
^  =  (cos(0  -I-  4>),  sm{e  +  (ji>)),  and  ^  =  (sin(0  +  (f>),  -  cos(0  -|-  (j>)). 


III.  SEPARATE  RECONSTRUCTION  OF  AND  Up(fc) 

In  this  section,  we  first  derive  a  frequency  domain  relationship  between  the  pro¬ 
jections  g(C,r)  and  scattering  potentials  C4(x)  and  Dp(x),  thus  obtaining  a  ”Pro- 
jection  Slice  Theorem”  (Deans,  1983)  associated  with  the  problem.  Inverting  this 
relationship  provides  both  a  frequency  domain  relationship  between  Up{}^, 

and  UB{k),  and  a  method  for  the  separate  reconstruction  of  ?7„(x)  and  Up{^,  or 
equivalently,  of  Ue{^  and  C/p(x). 

From  (27),  we  have 


g[k^,kr)  —  g^{k^,kr)  —  gp{k(,kr). 


(36) 
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where  K)  and  k^)  are  the  2-D  Fourier  transforms  of  g^ii^ »")  and  g^H,  r), 
respectively.  Since  =  g{^,r)  for  the  constant  density  case,  it  was  shown  in 

(Ozbek  and  Levy,  1987,  Section  4  and  Appendix  B)  that  g^{k^,kr)  can  be  expressed 
as 


L(*5.  M  =  -T=M  cr.  (k  =  kA+  +  7Sgn(/:,)^/*Flf  })  , 

V^r  -  ^  ^ 

(37) 

for  |A:f I  <  jkrl,  where 


^K(i)  =  /  dx'?7^(x')e-‘^-^' 


(38) 


is  the  2-D  Fourier  transform  of  and 


1=  \ 


-1-1  ifx-^— p>0  for  all  X  e  V , 
—1  if  X  •  ^  —  p  <  0  for  all  X  G  V . 


7  describes  which  side  of  the  receiver  array  the  support  V  of  the  inhomogeneities  is 
situated. 

We  now  express  p^(/:{,A:,.)  as  a  function  of  the  2-D  Fourier  transform  Up{l^  of 
Z7p(x^).  We  first  take  the  Fourier  transform  of  gp{^^r)  with  respect  to  r: 


/OO 

(irgpiC,  r)e“’‘'’' 

•OO 

=  f;u,k) 

=  (39) 

where  Fp  denotes  the  complex  conjugate  of  Fp  and  Ho^^(-)  is  the  Hankel  function 
of  order  0  and  type  2.  Now,  taking  the  Fourier  transform  with  respect  to  £  gives 


(40) 
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where  (Ozbek  and  Levy,  1987,  Appendix  B) 

N{^,k^,kr)  ^  r 

J  — OO  “* 

^  2sgn(<:,) 

for  |fcf|  <  \kr\.  Consequently,  we  get 


(41) 


ffp(^i,kr)-  =  - 


gi-IPsgnCk 


•t7p  (fc  =  krl  +  +  'tsgn(kr)^k^  -  k^  ^  , 

for  \k^\  <  |A:r|.  Combining  (36),  (37)  and  (42),  we  obtain 


+  '■isgn{kr)yjkl  -  A:| 


(42) 


kr)  = 

Up  (k  =  kri  +  +  Tfsgn(A;,.)y^A:2  -  A:|  ^ 


6  ■  +  •7sgn(fc,.)  ^k^  -  kl  $ 

kr 


(43) 


for  |A:f|  <  |A:,.|.  For  \k^\  >  |A;,.|,  g{k^,kr)  is  related  to  the  part  of  the  observed 
scattered  field  that  corresponds  to  evanescent  waves  (Ozbek  and  Levy,  1987),  and 
we  do  not  make  use  of  this  portion  of  g{k^,  kr)  in  our  inversion  scheme  .  The  inverse 
formula  of  (43)  is 


Ur{!^ 


Ucik)  -  2 


ik-iV 


k^ 


Up{k) 


2'Kk  •  6 


ipk-±/2k-l^g,  K  = 


A,  1 

k-t 

2k  ■  § 


,  kr  = 


k^ 


2k -9 


keC,  (44) 
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where  the  cone  C  is  defined  below,  and  where  Uc{k)  =  Uk.{1£)  +  U /,[]£)  is  the  2-D 
Fourier  transform  of  the  velocity  scattering  potential  ?/«(£)•  Here,  Ur{}^  denotes  the 
2-D  Fourier  transform  of  the  reconstructed  potential  Uji{^  obtained  by  applying 
the  constant  density  reconstruction  procedure  to  the  projections  r)  obtained 
from  a  variable  density  and  velocity  medium. 

Equation  (43)  represents  the  "Projection  Slice  Theorem”  associated  with  the 
variable  density  inverse  acoustic  problem  relating  the  1-D  Fourier  transform  of 
g{i,kr]  with  respect  to  ^  to  a  semicircular  slice  of  the  2-D  Fourier  transform  of 
Ur{^.  For  a  fixed  k,.,  g{k^,kr)  gives  UR{k)  along  a  semicircle  of  radius  |Ar,.|  centered 

A 

at  as  shown  in  Fig.  3.  By  letting  kr  vary,  these  semicircles  span  a  cone  C,  which 
is  defined  as 

C  =  {k:  \k'k\>  v/2/2}  (45) 

for  7  = +1,  where  =  (k,k)  and  ^  =  (cos((0  +  ^)/2),  sin((5  +  ^)/2)).  The  angular 
range  of  this  cone  is  90°,  as  indicated  in  Fig.  3.  For  7  =  —1,  C  is  the  complement 
C  of  the  above  cone. 

Combining  (34)  and  (44)  gives 

A 

UrUc)  =  Usik)  =  U^k)  -  2  cosV U,{k),  kGC.  (46) 

where  f  is  the  angle  between  the  vectors  fc  and  0,,  with  k=  (fc,  A).  This  relation  is 
the  key  result  of  our  paper.  It  shows  that  the  reconstructed  image  C)’ie(fc)  can  be 
obtained  by  applying  the  2-D  filter  'yX-rj}f2Tr^  to  the  backprojected  image  Usik)- 
This  relation  is  similar  to  the  identity  which  is  used  for  2-D  backprojection  and 
filtering  x-ray  reconstruction  methods  (Deans,  1983).  Since  the  filter  A  •  ^  diverges 
for  high  frequencies  k,  it  needs  to  be  "clipped”  as  k  becomes  large.  Furthermore, 
identity  (46)  shows  that  Ur{!£)  is  a  linear  combination  of  Uc{k)  and  Up{k).  This 
implies  that  it  is  not  possible  to  reconstruct  these  two  potentials  from  a  single 
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experiment.  To  reconstruct  Ue{k)  and  Up{!£i  separately,  we  need  in  principle  two 

A  A 

experiments  with,  plane  waves  incident  at  angles  ^  and  then  we  can  solve  the 
system 

1  -2{k-Ly 

MikLL) 

A  A  A 

which  requires  inverting  the  matrix  M(fc;^,£2)- 

For  the  numerical  stability  and  robustness  of  the  matrix  inversion  procedure, 
the  matrix  M{k^  6,1,0^)  must  be  as  nonsingular  as  possible.  The  most  appropriate 
measure  of  the  singularity  of  a  matrix  is  the  smallest  singular  value  of  the  matrix. 
The  smallest  singular  value  of  the  matrix  M  is 


UmiJc) 

UR2{k) 

(47) 


O^min  (■^^) 


min  a,-  such  that  detfcr?/  —  M'M)  —  0 

f=l,2  '  ' 


|277i  +  2772  +  1  - 


(277^  +  2772  +  1)^-4  (tjI  -  ntf 


1/2 


7  (48) 


^  A  A  \  ^  A 

where  771  =  it  •  £1  and  772  =  J£  -  0^. 

Inversion  of  M  would  be  most  robust  when  is  maximized.  This  takes 

AA  A  AA  AAAA 

place  for  values  of  kj0,i,  and  0^  such  that  0.1-02  =  0  and  =  ±£i  or  =  ±^. 
Therefore  the  two  probing  waves  must  be  incident  at  angles  perpendicular  to  each 
other.  Under  this  condition,  let  us  consider  the  frequency  domain  coverage  we 
would  have  for  finite  bandwidth  data,  assuming  that  we  have  receiver  coverage 
surrounding  the  medium.  Neglecting  the  low  frequency  cutoff  band,  the  frequency 
domain  coverage  due  to  a  single  probing  wave  has  a  "figure-of-eight”  shape  aligned 
with  the  direction  of  the  probing  wave  (Ozbek  and  Levy,  1987).  When  two  probing 
waves  are  used,  17(.(fc)  and  can  be  solved  only  in  regions  where  where  there  is 

double  coverage,  as  indicated  by  the  shaded  areas  in  Fig.  4.  However,  if  we  consider 
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the  superimposed  "radiation  pattern”  of  o'niin(M)  drawn  in  Fig.  4  also,  we  see  that 
M  is  most  singular  for  those  values  of  k  where  we  have  double  coverage. 

A  A 

For  general  values  of  i_i  and  0^  the  situation  is  similar.  In  general,  M  is  singular 

A  A  A  A  A  A  A 

for  values  of  which  satisfy  •  £.i|  =  |^  •  ^|.  Therefore,  for  =  ±^,  M  is  singular 
forall^  otherwise,  it  is  singular  for  fc  =  ±(^+^)/lli+^|  orfc  =  ±(£i— ^)/l£i— ^|. 
These  are  the  directions  which  in  fact  bisect  the  regions  where  there  is  double 
coverage.  The  coverage  and  radiation  pattern  for  the  case  when  the  angle  between 

A  A  ' 

and  0^  is  45°  is  shown  in  Fig.  5. 

The  preceding  analysis  confirms  what  researchers  in  diffraction  tomography  and 
exploration  geophysics  have  suspected  for  some  time:  namely,  that  it  is  exceedingly 
difficult  to  reconstruct  simultaneously  the  velocity  and  density  perturbations  of 
an  acotistic  medium  in  a  numerically  stable  way.  A  purely  theoretical  analysis 
neglecting  the  bandlimitation  of  the  probing  wave  would  lead  to  the  conclusion  that 
only  two  experiments  at  diflferent  angles  of  incidence  are  necessary.  Yet,  as  shown 
above,  the  numerical  results  obtained  by  such  an  approach  would  be  worthless.  In 
practice,  one  must  use  substantially  more  than  two  angles  of  incidence,  say  angles 

A  A  A 

and  for  each  k,  solve  the  resulting  system 


v'.W 

•  • 

* 

.1  -2{k-k,)\ 

-^(^)  djiik) 


(49) 


by  the  least  squares  method,  where  {ty,  zy , . . . ,  c  {1,2,  ...,iV}  is  the  set  of 

indices  corresponding  to  the  angles  of  incidence  for  which  the  probing  wave  provides 
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coverage  at  k_.  This  gives  the  solution 


Ucik) 


IV.  THREE-DIMENSIONAL  GEOMETRY 


After  discussing  the  2-D  experimental  geometry,  we  summarize  the  correspond¬ 
ing  results  for  the  3-D  case.  For  the  3-D  geometry,  we  assume  that  the  receivers  are 
on  a  plane;  for  convenience  we  choose  this  to  be  the  x-y  plane.  The  position  of  an 
arbitrary  receiver  located  in  this  plane  is  therefore  given  by  £  =  (^,0),  where  ^ 
represents  the  x-y  coordinates  of  the  receiver.  The  Green’s  function  due  to  a  point 


source  is 


Goix,3^,u>)  = 


^ik\x-x‘\ 


4ir\x-^\' 

Under  the  Born  approximation,  the  Lippmann-Schwinger  equation  takes  the  form 
(19),  and  the  projections  in  this  case  become 

=  -  ff^{|y,r),  (52) 


where 


fdx'U(T')\ff  ,  sgn(r-|.i'-|x'- ID 

/  -  "'-U-  li^J  i  iiPa  J’ 
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where  0  is  the  unit  vector  in  the  direction  of  propagation  of  the  incident  wave.  The 
projection  r)  is  again  a  surface  integral  of  the  scattering  potential  over  a  circu¬ 

lar  paraboloid,  like  in  the  constant  density  case.  The  weighting  function  appearing 
in  the  representation  (53)  of  is  an  impulse,  and  in  this  sense  the  3-D  case 

is  quite  different  from  the  2-D  case.  Projection  gp{^^,r)  shows  further  deviations. 
Like  in  the  2-D  case,  it  has  a  cosine  dependence  on  the  angle  between  the  direction 
of  the  incident  wave  and  the  direction  of  scattering.  The  weighting  function,  in 
addition  to  an  impulse,  contains  a  smoother  term  which  in  fact  makes  gp{^^,r)  a 
noncausal  function  of  r.  This  is  due  to  the  l/k^  filter  that  is  applied  to  the  scat¬ 
tered  field  P,(£,  w)  to  obtain  the  projection  ffp(|^,r).  If  the  far-field/high-frequency 
approximation  A:|x'  —  1  is  made,  the  impulsive  term  clearly  dominates. 

We  introduce  the  backprojection  operation  as  in  the  constant  density  case: 


=  J  j^^drg{^^,r) 


5(r  —  £•  X  —  |x  —  £|) 

k-  ^1 


(55) 


In  the  frequency  domain,  the  backprojected  image  can  be  related  to  the  parabolical 
projections  as  (see  Ozbek  and  Levy,  1987,  Section  8  and  Appendix  C) 


Ueik)  = 


12%  A 


A.-I  X.-i\ 
2k-l’2k-ij 


(56) 


where 

kktT,  kr)  =  I  drg((^,  r)e-‘(i<r 

are  the  3-D  Fourier  transforms  of  Ub{x}  and  ff(^y,,r),  and 


A,  4  {kl-kl-  kl  2k,k„  2k^k.), 
A,  4  (2k,k„  kl-kl -kl2k,k.). 


(57) 

(58) 
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From  (52),  we  again  have 

A  A  A 

dik^T^^r)  —  “  9p{k^Ty  ^r)y 


(59) 


A  A 

where  and  5p(^T>^r)  are  the  3-D  Fourier  transforms  of  r)  and 

9p{ij.,r),  respectively.  Again,  since  3«(|y,r)  =  9{^,r)  for  the  constant  density 
case,  it  was  shown  in  (Ozbek  and  Levy,  1987,  Section  8  and  Appendix  D)  that 

A. 

takes  the  form 

hik^T,  K)  =  -  (-  ^  -  I^tH)  ,  (60) 

for  l^xl  ^  l^rl)  where  U^ik)  is  the  3-D  Fourier  transform  of  U^{^,  and 


7  =  < 


+  1  if  2  >  0  for  all  X  6  V, 
—1  if  2;  <  0  for  all  X  G  V. 


Using  the  intermediate  results  derived  in  (Ozbek  and  Levy,  1987,  Appendix  D), 
we  also  have 


hik^T,kr)  =  - 


i2jr 


i-  {^r>7sgn(A:r)\/A:2  -  |i^Tp)] 


\^r\\Jkl  —  l^rl 
•Up  (fc  =  {kiT,15gn{kr)\Jkl  -  l^rp))  , 

for  l^rl  <  |A:,.|,  where  Up{l^  is  the  3-D  Fourier  transform  of  Up{^. 
Combining  (52),  (59)  and  (60)  yields 


(61) 


g{k^Tikr)  — 
i2irsgn{kr) 

\Jkl  —  l^rl' 


{Uk  (k  =  kri+  (^r,7sgn(fcr)\/fe?^HSrF)) 


i-  {k^T,lsgn(k,)yJkf  -  |^j,|2)  - 

-  Up  (fc  =  kri  +  (^r.  l^gri{kr)\Jkl  -  }  , 

(62) 
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for 

The  inverse  formula  of  (61)  is 


UrUc) 


A 


UciJc)  -  2 

•  i  i. 


ik-ir 


uAk) 

kii  vi) 

2k -i'  2k -ij 


kr 


k^C, 

(63) 


where 


2k^k^,  kl-kl-  kl), 


and  where  UdM)  =  U^k)  +  Up{k)  is  the  3-D  Fourier  transform  of  the  velocity- 
scattering  potential  Ue{x).  Here,  Ur{Jc1  denotes  the  3-D  Fourier  transform  of  the 
reconstructed  potential  Ur{x)  which  is  obtained  by  applying  the  constant  density 
reconstruction  procedure  to  the  projections  g{^,  r)  obtained  from  a  variable  density 
and  velocity  medium. 

For  l^rl  >  i^r|>  S'S  in  the  2-D  case,  g(^r>^r)  is  related  to  the  part  of  the 
observed  scattered  field  that  corresponds  to  evanescent  waves,  and  we  do  not  make 
use  of  this  portion  of  g{k^Ti  kr)  in  our  inversion  scheme. 

Equation  (61)  represents  the  "Projection  Slice  Theorem”  for  the  3-D  case  as¬ 
sociated  with  the  variable  density  inverse  acoustic  problem.  For  a  fixed  kr,  the 
2-D  Fourier  transform  of  g{^j.,kr)  gives  the  3-D  Fourier  transform  of  17r(x)  over  a 
hemisphere  of  radius  |A:r|  centered  at  kri..  By  letting  kr  vary,  U{k)  is  determined  in 
a  cone  C,  which  again  covers  half  of  the  3-D  frequency  space. 

Combining  (56)  and  (62)  gives 


=  =  tec. 


(64) 
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Therefore,  Ur{^  for  JceC  can  be  obtained  from  Usik)  by  a  3-D  space  invariant 
filtering  operation. 

Ur{I^  is  related  to  Udk)  and  Up{!c)  with  the  same  relationship  ((46)  and  (64)) 
in  two  and  three  dimensions.  Therefore,  the  procedure  for  individual  reconstruction 
of  ?7c(x)  and  £/p(x)  is  the  same  in  these  two  cases. 

V.  NUMERICAL  EXAMPLE 

The  theory  presented  in  this  paper  was  tested  for  the  two-dimensional  case,  using 
computer-generated  synthetic  data.  Figs.  6  and  7  show  the  velocity  and  density 
scattering  potential  models,  Uc{x)  and  17p(x),  respectively.  The  scattering  potentials 
correspond  to  velocity  and  density  anomalies  which  are  constant  in  square-shaped 
areas  of  dimensions  35  m  x  35  m.  The  background  medium  was  homogeneous 
with  a  velocity  of  5000  m/s  and  a  density  of  2000  kg/m^.  The  magnitude  of  the 
velocity  and  density  perturbations  were  10%  of  the  nominal  values.  A  lowpass 
source  wavelet  with  a  cutoff  frequency  of  425  Hz  was  used,  so  that  the  object  sizes 
are  three  times  the  shortest  wavelength  in  the  source  signal.  The  regions  of  anomaly 
are  separated  by  a  distance  equal  to  six  times  the  shortest  wavelength.  The  synthetic 
scattered  waves  were  obtained  by  using  the  forward  scattering  equation  under  the 
Born  approximation;  however,  since  the  velocity  and  density  perturbations  are  of 
limited  size  with  respect  to  the  shortest  wavelength,  and  the  magnitudes  of  these 
perturbations  are  relatively  small,  we  do  not  feel  that  the  approximation  is  critical 
for  this  particular  example.  The  entire  image  area  was  500  m  x  500  m,  the  grid 
size  was  5  m  x  5  m,  and  receivers  were  located  on  all  sides  around  the  medium, 
100  on  each  side. 

As  indicated  in  Section  III  of  this  paper,  in  order  to  guarantee  the  numerical 
stability  of  the  procedure  for  reconstructing  separately  the  velocity  and  density’ 
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inhomogeneities,  more  than  two  sources  are  needed.  In  this  experiment,  we  have 
used  eight  angles  of  incidence,  at  22.5®  intervals.  The  inversion  was  performed  over 
the  regions  in  the  k  domain  where  coverage  was  provided  by  at  least  five  probing 
waves;  i.e.,  using  the  notation  of  eq.  (49),  N  —  8,  P  >  5,  and  rank(M)  >  4  for  all 
inversion  points  k.  This  corresponds  to  carrying  out  the  inversion  over  a  circular 
lowpass  region  with  a  radius  of  about  55%  of  the  maximum  frequency  coverage 
provided  by  a  single  source. 

For  comparison,  we  first  examine  images  and  U^{^.  To  obtain  these 

images,  in  the  frequency  domain,  values  obtained  due  to  different  sources  providing 
multiple  coverage  were  simply  averaged  point  by  point.  Fig.  8  shows  the  back- 
projected  image  17^" (x)  can  be  interpreted  as  a  "migrated”  image  of  the 

velocity  field  for  a  constant  density  medium  (for  migration,  see,  for  example,  Claer- 
bout,  1985).  Fig.  9  depicts  Z7^"(x),  which  is  the  image  obtained  by  applying  the 
constant  density  reconstruction  procedure  to  the  data  obtained  from  a  variable 
density  and  velocity  medium. 

Some  observations  can  be  made  regarding  these  images.  Both  images  display  the 
locations  of  the  scatterers;  however  the  "inversion”  image  is  much  better  focused 
than  the  ’’migration”  image.  This  effect  has  also  been  noted  by  other  researchers 
(Esmersoy  and  Miller,  1987).  In  addition,  the  values  of  I7^”(x)  differ  by  orders  of 
magnitude  from  the  numerical  values  of  the  true  scattering  potentials.  On  the  other 
hand,  looks  like  Ue{s^  —  Up{g^,  and  actual  constructed  values  confirm  this.  To 

interpret  this  result,  observe  from  (46)  that,  for  the  "ideal  case”  where  plane  wave 
experiments  are  performed  for  all  angles  $  of  incidence,  the  averaging  scheme  de¬ 
scribed  above  for  combining  the  reconstructed  images  obtained  for  different  probing 
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angles  can  be  written  as 

Uk"  m  =  l  r  MUn(.k)  =  p,(i)  -  U,m  =  Cf.ca,  (65) 

where  Ur  corresponds  to  the  inversion  result  for  one  source,  and  corresponds 
to  the  result  after  angular  averaging.  Therefore,  averaging  the  reconstructed  poten¬ 
tial  Ur{I^  over  different  angles  is  equivalent  to  reconstructing  the  compressibility 
potential  of  the  medium. 

Figs.  10  and  11  show  the  separate  reconstructions  of  Uc{^  and  Up[^,  respec¬ 
tively.  The  numerical  values  obtained  are  within  20%  of  the  model  values.  There 
exist  several  sources  of  error.  The  first  of  these  is  the  bilinear  interpolation  proce¬ 
dure  which  is  used  to  convert  the  discretized  2-D  Fourier  transform  g{k^,  of  the 
projections  into  the  discretized  Fourier  transform  Ur{1z)  of  the  reconstructed  image. 
A  second  source  of  error  is  the  fact  that  all  the  inversion  results  developed  in  this 
paper  assume  that  receiver  arrays  are  infinite,  whereas  the  arrays  which  are  used 
for  the  present  example  have  a  finite  length.  Other  errors  are  due  to  the  fact  that 
the  source  wavelet  is  bandlimited,  and  needs  to  be  deconvolved. 

Finally,  the  DC  levels  of  the  velocity  and  density  scattering  potentials  cannot  be 
reconstructed  separately  with  the  derived  inversion  formulas,  since  the  coefficient 
of  Up{]^  in  equation  (44)  is  not  analytic  around  k  =  0.  The  DC  level  of  can 

be  recovered  from  equation  (22),  so  that  if  the  DC  level  of  either  the  density  or  the 
velocity  is  known,  the  other  one  can  be  computed.  In  our  implementation,  we  have 
estimated  (fc  =  0)  and  Up^k  —  Oi)  as  a  weighted  average  of  the  closest  eight  values 
in  the  discrete  wavenumber  domain.  Adopting  the  reciprocal  of  the  square  of  the 
distance  as  the  measure  of  weight,  we  assigned  a  weight  of  1/6  to  the  closest  four 
samples,  and  1/12  to  the  next  closest  samples  which  are  diagonally  located. 
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VI.  CONCLUSIONS 

In  this  paper  we  have  considered  the  problem  of  the  reconstructing  separately 
of  the  velocity  and  density  inhomogeneities  of  a  multidimensional  acoustic  medium 
probed  by  wide-band  plane  waves.  The  problem  was  posed  as  a  generalized  tomo¬ 
graphic  problem,  where  weighted  integrals  of  the  velocity  and  density  scattering 
potentials  Uc(x)  and  Up(x)  are  used  as  data.  A  backprojection  operator  C/s{z)  was 
defined,  which  was  related  to  the  generalized  projections  in  the  Fourier  transform 
domain.  It  was  shown  that,  by  applying  a  time-invariant  filter  to  Z7b(x},  we  can 
obtain  an  image,  Ur{^,  which  in  the  Fourier  domain  is  a  linear  combination  of 
the  velocity  and  density  scattering  potentials,  and  where  the  coefiicients  depend 
on  the  angle  of  incidence  of  the  probing  wave.  Therefore,  for  numerical  stability, 
several  angles  of  incidence  were  used  to  solve  for  the  velocity  and  density  scattering 
potentials  separately. 
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Figure  Captions 

Fig.  1  2-D  experimental  geometry. 

Fig.  2a  Scattering  pattern  due  to  velocity  inhomogeneities. 

Fig.  2b  Scattering  pattern  due  to  density  inhomogeneities. 

Fig.  3  Coverage  of  Ur{^  for  a  single  array. 

A  A  A 

Fig.  4  Frequency  coverage  of  Ur{^  and  the  "radiation  pattern”  of  <Tnun(M;  k;  £i,  £2) 
for  the  case  when  two  probing  waves  are  used,  incident  at  right  angles  to  each 
other. 

Fig.  5  Frequency  coverage  when  of  Ur{J^  and  the  "radiation  pattern”  of  cr„un{^',  £i)  £2) 

for  the  case  when  two  probing  waves  are  xised,  incident  at  a  45°  angle  with 
respect  to  each  other. 

Fig.  6  Velocity  scattering  potential  model  for  the  synthetic  experiment. 

Fig.  7  Density  scattering  potential  model  for  the  synthetic  experiment. 

Fig.  8  The  backp rejected  image  Ur’{x),  obtained  assuming  a  constant  density 
medium. 

Fig.  9  The  inverted  image  obtained  assuming  a  constant  density  medium. 

Fig.  10  The  separate  reconstruction  of  the  velocity  scattering  potential. 

Fig.  11  The  separate  reconstruction  of  the  density  scattering  potential. 
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